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Aims. We present a new C++ code for collisional N-body simulations of star clusters. 

Methods. The code uses the Hermite fourth-order scheme with block time steps, to advance the particles in time, 
while the forces and neighboring particles are computed using the GRAPE-6 board. Special treatment is used for close 
encounters and binary or multiple subsystems that form either dynamically or exist in the initial configuration. The 
structure of the code is modular and allows the appropriate treatment of more physical phenomena, such as stellar and 
binary evolution, stellar collisions, and evolution of close black-hole binaries. Moreover, it can be easily modified so that 
the part of the code that uses GRAPE-6, could be replaced by another module that uses other accelerating-hardware 
such as the graphics processing units (GPUs). 

Results. Appropriate choice of the free parameters give a good accuracy and speed for simulations of star clusters up 
to and beyond core collapse. Simulations of Plummer models consisting of equal-mass stars reached core collapse at 
t :^ 17 half-mass relaxation times, which compares very well with existing results, while the cumulative relative error 
in the energy remained below 10~^. Comparisons with published results of other codes for the time of core collapse for 
different initial conditions, show excellent agreement. Simulations of King models with an initial mass-function, similar 
to those found in the literature, reached core collapse at t :^ 0.17 ± 0.05 half-mass relaxation times, which is slightly 
earlier than expected from previous work. Finally, the code accuracy becomes comparable to and even better than 
the accuracy of existing codes, when a number of close binary systems is dynamically created in a simulation. This is 
because of the high accuracy of the method that is used for close binary and multiple subsystems. 
Conclusions. The code can be used for evolving star clusters containing equal-mass stars or star clusters with an initial 
mass function (IMF) containing an intermediate mass black hole (IMBH) at the center and/or a fraction of primordial 
binaries, which are systems of particular astrophysical interest. 
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1. Introduction 

Systems of numerous individual members interacting with 
each other are found on both small and large scales. Gas 
molecules interact with each other by means of electrostatic 
van der Waals forces, the same forces are responsible for the 
interactions and collisions of amino-acids that lead to the 
construction of protein molecules. On the other hand, plan- 
ets and planetesimals interact with both each other and the 
central star in a planetary system, with the gravitational 
force. The gravitational force is also the interaction force 
between members of greater systems such as open or glob- 
ular star-clusters, galaxies, and groups of galaxies. 

N-body codes provide one way of simulating systems 
like those mentioned above. Those codes require calcula- 
tions of all the forces between all the individual particles 
at every integration time step. This means that at every 
time-step, a total of N{N — 1) ^ N"^ force calculations 
should be made, where N is the number of interacting par- 
ticles. From that, it is easy to conclude that N-body codes 
have time complexity 0{N'^) and require fast computers or 
even special purpose hardware, such as GRAPE systems, 
to simulate realistically large systems. This is why the evo- 
lution of N-body codes that could integrate these systems 
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in a reasonable time, is closely connected to the hardware 
evolution. 

The first attempt to perform a simulation, albeit with- 
out the use of a computer, but instead a system of 37 
light bulbs represending two interacting model galaxies, 
was made by Holmberg Holmberg (1941). The first com- 
puter N-body simulations were developed by von Hoerner 
(1960, 1963) using initially 16 and later 25 particles. Later, 
as simulation methods and computers became faster, the 
number of particles increased to 100 (Aarseth 1963) and 
250 (Aarseth 1968). Even in these early attempts some of 
the physical processes that occur in true star clusters, such 
as mass segregation and formation of hard binaries, were 
found in the simulations, but simulating systems closer to 
realistic star clusters required improvements in both soft- 
ware, by developing efficient algorithms, and hardware. 

A significant step in the evolution of N-body algorithms 
was the introduction of individual time steps: all parti- 
cles do not share the same time step, but each one of 
them has its own, depending on its motion. In block-time- 
step schemes (McMillan 1985), the time steps are quan- 
tized, usually as powers of 2, permitting blocks of par- 
ticles to be updated simultaneously. This idea made the 
algorithms faster, because only the necessary calculations 
need to be made for a given accuracy. Today all sophis- 
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ticated N-body codes such as Starlab^ (Portegies Zwart 
et al. 2001; Hut 2003), ^GRAPE (Harfst et al. 2008), 
and Aarseth's nbody^ series (nbody4, nbody6, nbody6++ 
(Aarseth 1999a, b; Spurzem 1999)) use block time steps. 
Another niilestone in the history of N-body simulations 
was the introduction of special purpose hardware, responsi- 
ble only for the calculation of the gravitational forces. The 
GRAPE series of computers (Makino et al. 1997, 2003a) 
were designed to perform only this type of calculation. The 
next step was GRAPE-DR^ (Makino 2008), which IS a more 
general purpose hardware specialized for N-body codes 
which operates faster than a normal CPU. Simulations of 
realistically large star systems can be made using a single 
GRAPE machine attached to a fast host-computer. Today, 
GPUs are slowly replacing CPUs and the older GRAPEs in 
N-body simulations (Portegies Zwart et al. 2007a; Belleman 
et al. 2008; Hamada & litaka 2007), because they are faster 
and considerably cheaper, than either CPUs or GRAPEs. 

In this paper, we assess the limitations of N-body codes 
with today's software and hardware power. First of all, N- 
body simulations are eithercollisionless or collisional. In col- 
lisionless simulations, the gravitational law is modified in 
close encounters to avoid collisions of particles and disen- 
courage the dynamical formation of hard binary systems. 
The modified gravitational force acting on a single particle 
i due to all other particles is calculated to be: 
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where e is the so-called softening parameter. The presence 
of close binary systems slows down the simulations, thus 
when handling large numbers of particles, one prefers to 
use collisionless simulations. 

In simulations of collisional systems, the minimum num- 
ber of stars that can be modeled is much lower and depends 
strongly on the fraction of stars that belong to binary sys- 
tems in the initial configuration (primordial binaries) . This 
is because the evolution of binary systems requires special 
treatment, including small time steps. Baumgardt et al. 
(2005) presented the results of simulations of unequal-mass 
star clusters containing up to A/" = 131072 stars and an 
intermediate-mass black hole (IMBH) at the center. Those 
simulations were oparated for 12 Gyrs and no primordial 
binaries were included. In Portegies Zwart et al. (2007b), 
the results of simulations of star clusters with a total of 

= 144179 stars including 10% primordial binaries were 
presented. The simulations were run for 100 Myrs. From 
the above numbers, it is easy to conclude that although it 
is not possible to simulate a typical galaxy with all its stars 
(AT ?^ 10^ — 10^), a simulation of a typical globular cluster 
{N ^ 10^ — 10^) for timescales of the order of some Gyrs is 
possible, with today's computing power. If primordial bina- 
ries are not included, the evolution of a star cluster with up 
to 10^ stars up to and beyond core collapse is possible, using 
a GRAPE-6 machine attached to a fast host-computer. In 
the near future, using new algorithms and the coming gen- 
erations of hardware, it should become possible to simulate 
the dense and massive centers of galaxies {N ^ 10^ — lO''). 



^ http://www.ids.ias.edu/~starlab/ also see http:// 
muse . li for a set of N-body simulation tools 

^ http : / / www . ast . cam . ac . uk/~ sverre/ web/ pages/nbody . 
htm 

^ http://www.kf cr.jp/index- e .html 



During the past few years, N-body simulations have 
provided a realistic tool for explaining the dynamics and 
structure of star clusters. Portegies Zwart et al. (1999) and 
Portegies Zwart & McMillan (2003) explained one way of 
forming IMBH in some globular clusters, by means of the 
process of runaway collisions of stars. Baumgardt et al. 
(2005) used detailed N-body simulations to discuss the 
structure of a globular cluster containing an IMBH, while 
Gualandris & Merritt (2008) discussed the possibility of 
the ejection of supermassive black holes from galactic cen- 
ters. Kupi et al. (2006) and Amaro-Seoane et al. (2009) 
used post-Newtonian effects on their N-body simulations 
to study the dynamics and the gravitational radiation from 
IMBH-binaries in star clusters. Finally, Baumgardt et al. 
(2008) discussed the effects of primordial mass segregation 
on the evolution of globular clusters. 

Here we introduce Myriad^, a new C++ N-body code 
for collisional and collisionless simulations. The code is 
loosely based on the instructive on-line books found in The 
Art of Computational Science website^. It has many of the 
features of existing codes, such as the use of block time 
steps, the Hermite 4*^ order integrator, GRAPE-6 support, 
and special treatment for binary and multiple sub-systems, 
but introduces a new structure and some new ideas espe- 
cially in using GRAPE-6 to find neighboring stars and per- 
turbers of binary systems. 

These ideas and a description of the basic structure of 
the code are presented in Sect. 2. In Sect. 2.2, the detailed 
steps of the code are described. More attention is given in 
describing the method for dynamical evolution of binary or 
multiple sub-systems that may appear in a simulation, in 
Sect. 2.2.3. In Sect. 3 we compare our preliminary results 
with the results of other codes. In the final two sections, 
some of the future directions of this code are presented with 
a discussion about the source of numerical errors and the 
CPU time needed for the simulations. The code uses the 
so-called N-body units, which, together with some basic 
equations for calculating the physical parameters of a star 
cluster are given in Appendix A. 



2. The new code 

2.1. Code structure 

The core of Myriad^ uses C++ classes to store and ma- 
nipulate data. The simplest or lower-level class is the class 
PARTICLE, which contains all necessary information about 
a single star^ (mass, radius, position, velocity, acceleration 
etc). The class particle also contains a set of functions 
that act on a star following its evolution in time, finding 
its neighbors, its time step and all the required data for 
the star itself. A set of particle classes that share the 
same time step at a given time, forms a higher-level class, 
the class BLOCK. As with the class particle, for the class 



^ Myriad is an ancient Greek unit for 10 000, which is the order 
of magnitude of the size of a star cluster that can be evolved 
using the code, in reasonable time. In modern English, the word 
refers to an unspecified large number http://eii.wikipedia. 
org/wiki/Myriad 

^ http://www.artcompsci.org 

^ Myriad will be soon available for download on http : //www . 
astro . auth . gr/~simos/mediawiki- 1.6. 7/index . php/Myriad 

^ A star is also referred as particle 
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CLUSTER 



• Total Mass 

• Half Mass Radius 

• Core Radius 

• Core Mass 

• Half Mass Relaxation Time 



BLOCK 



• Time Step 

• Number of Stars 
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• Mass 

• Time Step 

• Position 

• Velocity 
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• Higher Derivatives 



Fig. 1. Class hierarchy in Myriad. A class CLUSTER is a set of several block classes. A block contains a variable number of 
PARTICLE and BINARY classes. A binary is a set of two or more particle classes. A class particle contains information (mass, 
position, velocity etc) for a single star and functions that act on them. The class cluster contains information about the system 
as a whole (total mass, energy, half-mass radius, core radius etc) and the required functions for calculating them. 



Starlab 



Initial Data 




Snai 


jshot 


Oul 


put 



Fig. 2. Simplified graphical representation of Myriad. The arrows show the flow of the data in the code. Boxes represent input or 
output datafiles or outside applications, while circles represent sets of code- functions. Myriad applications are those inside the large 
dashed box. The other boxes and circles refer to satellite programs that create initial data or make animations from Myriad-output. 



BLOCK there is a set of functions that act on both the class 
itself and its particle members. 

Finally, the complete set of BLOCK classes (it might only 
be one if all stars evolve with the same time step) forms a 
cluster, which is the highest-level class. The CLUSTER- 
functions, act on all particle classes of the system, and 
mainly finding information about the cluster of stars itself 
(total mass, center of mass, escapers, binaries etc). There 
is also another class that is a member of a class BLOCK, 
the class binary, which consists of two or more parti- 



cle classes representing stars that lie close to each other in 
space. The binary class has the required functions for the 
initialization of a binary or multiple star-system, its evolu- 
tion in time, and its termination, when necessary, and also 
returns the required data back to the CLUSTER. 

Figure 1 shows the hierarchy of classes in the code de- 
scribed above. Figure 2 is a simplified graphical represen- 
tation of Myriad. Initial data are provided to Myriad from 
a single data file. This initial data files can be constructed 
by the Starlab package. The main output of Myriad are 
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snapshots of the evolving system at every output time step 
dt out and information about the whole cluster (total mass, 
core radius, half-mass radius etc) at every diagnostic time 
step dt diag. 

2.2. Integration method 

The integration method used in Myriad is the Hermite 4^^ 
order (H4) (Makino & Aarseth 1992) predictor-corrector 
scheme (PEC: Predict-Evaluate force-Correct) with block 
time steps. Accelerations and their derivatives (jerks) are 
computed using GRAPE-6. Binaries, close encounters, and 
multiple sub-systems that require small time steps are 
detected using GRAPE-6 and evolved using the time- 
symmetric Hermite 4*^ order integrator (Kokubo et al. 
1998), where the prediction (P) step is the same as in the 
simple H4, but the force evaluation (E) and correction (C) 
steps are repeated 3 times {P{EC)^). 

For each BLOCK i, ti^ is the current time of its particle 
classes, dti is the time step, and t^ = tic + dti is the time 
after a single update. 

The integration procedure consists of the following steps 
(the time complexity of the slowest steps is displayed in 
parenthesis) : 

1. We define the universal current time tc and the time 
after the step is completed tf to be equal to the time 
and time forward of BLOCK 0, the BLOCK containing 
stars with the minimum time step, respectively 



(c) We evaluate the accelerations and their derivatives 
for all the particles of the BLOCK using the GRAPE- 

6 hardware 

(d) We find the new time step for each particle of 
the BLOCK. This will be used for the next update 

(0(m)). 

(e) We check for encounters between single stars or bi- 
naries. 

(f) We correct the positions and velocities of all the par- 
ticles of the BLOCK (^0{m)^. 

(g) If an encounter is detected in step (e), we construct 
the BINARY, and replace its members with a "vir- 
tual" star represending their center of mass. 

(h) We check the binary classes. If a binary is termi- 
nated, replace its center of mass in the N-body code, 
with its former members and erase it. 

(i) We then send the corrected values to GRAPE-6 
memory. 

4. We move particles between BLOCK classes according to 
their new time step. 

5. We update tc and tf for each updated BLOCK. 

6. We continue with step 1. 

From the above, it is obvious that the speed of the code 
depends on: 



tc — toc, (2) 

tf = tof. (3) 

2. We determine which BLOCKS have time forward equal 
to tf and label them as requiring update. 

3. We update each labeled BLOCK to its time forward. 
To update a single BLOCK with m members (A^ being 
the total number of stars): 

(a) If binary^ classes exist, update them to tf 
To update a binary class with k members: 

i. We determine the time step dt^ for the P{EC)^ 
integrator. 

ii. We predict the position and velocity of each 
member (^0{k)^. 

iii. We calculate the forces between binary- 
members on the CPU (^0{k'^)^. 

iv. We calculate the perturbation forces from p 
neighboring stars on the CPU (^0{kp)^ . 

V. We then correct the position and velocity of each 
member (^0{k)^. 

Processes (iii) to (v) are repeated n = 3 times. 

vi. We update the current time tc = + dt^ of the 

BINARY. 

vii. We continue from (i) until tc equals tf. 

(b) We predict the positions and velocities of all the 
particles of the BLOCK, including normal stars and 
"virtual" stars that represent the centers of mass of 

binary or multiple sub-systems ^O(m)^ . 
^ A BINARY can contain two or more stars. 



1. The total number of stars N. 

2. The number of BLOCK classes and the average number 
of stars per BLOCK. 

3. The average number p of stars responsible for significant 
perturbations on a binary system (perturbers) . 



Figure 3 shows how data flows between the CPU and 
the GRAPE-6 and which parts of the calculation are per- 
formed on the CPU and the GRAPE-6. Initially positions 
and velocities, and estimates of both the accelerations and 
the jerks of all particles are stored in GRAPE-6 memory 
(This procedure is omitted in the graph). Then, a predic- 
tion is made for a set of i-particles using the CPU. The 
predicted values and the time of each particle is sent to the 
GRAPE-6 buffer, which directly sends them to the inter- 
nal predictor of the GRAPE-6. The internal predictor takes 
the data for the rest j-particles from the GRAPE-6 memory 
and predicts their positions and velocities at the time of the 
i-particles. The GRAPE-6 then computes the acceleration, 
its derivative (jerk), and the nearest neighbors for the i- 
particles and returns them to the corrector operating on the 
CPU. The corrector calculates the corrected values of posi- 
tion and velocity for the i-particles and sends the new data 
to the GRAPE-6 memory. A new set of i-particles then, is 
sent to the predictor and the procedure continues until all 
particles are evolved to the expected time. Data flow be- 
tween the CPU and the GRAPE-6 occurs three times per 
time step. The most time-consuming calculation, the calcu- 
lation of the gravitational forces {0{N log{N)))^ performed 
on the GRAPE-6. The rest of the calculations are of order 
0{N). If the force-calculation were to be performed on the 
CPU, then it would be far more time-expensive as it would 
scale as 0{N'^). 
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Fig. 3. Schematic data flow between the CPU and the GRAPE-6 for the Hermite 4*^-order integrator. See text for description. 
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Fig. 4. Distribution of equal mass particles in different time steps (powers of 2). Aarseth's criterion is used for the time step with 
r] — 0.01. A small number of particles lies in time step 2~^^. This is the smallest allowed time step for the N-body code. Stars 
with this time step are candidates for close encounters or virtual particles representing the centers of mass of binary or multiple 
sub-systems 



2.2.1. Block time steps 

The use of block time steps have proven to be an ideal 
method for N-body simulations. The advantage of block 
time steps over constant (shared) time steps for all particles 
is that a smaller number of operations is needed to achieve 
a given accuracy. The advantage over individual time steps 
is that particles are monitored with time in groups and 
not individually. The implementation of block steps is ideal 
with GRAPE-6, because it lessens the communication time 
between the host and GRAPE-6. In individual time-step 
algorithms, a single particle is sent to GRAPE-6 every time, 
while in block time step algorithms, groups of particles. 



that share the same time step are sent simultaneously to 
GRAPE-6 for force evaluation. 

At t = 0, all particles, unless there are binaries in the 
initial configuration, share the same time step, which is the 
smallest time step allowed by the N-body code. This time 
step is given by (Aarseth 2003) 

Dt^.'OM{^) (i) , (4) 

where fa is the mean mass of the system, 77/ the initial 
accuracy parameter with typical value 0.01, and R^i is the 
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close encounter distance 



2.2.2. The Hermite 4*^ order predictor-corrector scheme 



^cl = 



2Gm 



(5) 



where G is the gravitational constant, a is the rms velocity 
dispersion given by 



GNm 



(6) 



at equilibrium, TV is the number of stars, and Ry is the 
virial radius of the cluster 



Ry = G 



N^fh 



(7) 



where U is the total potential energy. 

The time step Dtmin is rounded to the closest power of 
2 and the resulting Atmin is the time step of BLOCK and 
the particles that belong to it. 

Each time a particle i is updated, its next time step is 
calculated according to Aarseth's criterion (Aarseth 2003) 



dtiA 



|ai,i||ai,i| + |ai,i| 



(8) 



where rj is the accuracy parameter with a typical value 0.01, 
and ai^i, ai^i, ai^i, and a'i^i are the acceleration of parti- 
cle i and its time derivatives, calculated in next section. 
Subscript "1" refers to the values of the parameters at the 
end of the time step (final values), while subscript "0" refers 
to the values at the beginning (initial values). 

The second and third time derivatives of the accelera- 
tion at the end of the time step are given by 



and 



(9) 
(10) 



where At^^o is the previous time step of particle i. 

Time steps are quantized according to the rule At^^i = 
2^, where n is a negative integer defined in such a way that 



2"" <dtii <2^+^ 



(11) 



In this way, all particles have time steps that are powers of 
2. All particles with the same time step are grouped in a 
BLOCK so that they are updated simultaneously. The time 
step of "virtual" particles replacing binaries or multiple sys- 
tems (see below) is set to be Atmin- In addition, Atmin is 
the time step of particles that are relatively close to and 
approaching each other rapidly. To improve the efficiency 
of the code, every new time step cannot be longer than 2 
times its previous value, jumps to higher time steps with 
step greater than 2 are not allowed. 

The diagram of Figure 4 shows an example of a particle 
distribution in different time steps in a cluster of 16 384 
(2^^^) stars. 



The H4 scheme is the integrator used in many N-body 
codes. The most important feature of this integrator is 
that it needs only one calculation of the acceleration and 
its derivatives (force calculation) in every time step. This 
makes H4 superior to other popular integrators of the same 
accuracy, which require more than one force calculations 
per time step. A Runge-Kutta 4*^ order algorithm evalu- 
ates the force three times every time step, which makes 
this integrator much slower than H4 and consequently in- 
appropriate for N-body simulations with large number of 
particles. 

Another reason for choosing the H4 integrator is that 
GRAPE-6 is specifically designed for use with this kind of 
integrator. This is because GRAPE-6 returns the accelera- 
tion and its first derivative for each particle, which are the 
only data H4 needs. One may use a higher order integra- 
tor of the Hermite family (Nitadori & Makino 2008), but 
this would require extra calculations to be made of higher 
order derivatives of the acceleration on the CPU. Higher 
accuracy, also is not necessary for this kind of problem. 

The integration steps for advancing a particle forward 
in time with the H4 scheme are given in the following lines 
(see also Makino & Aarseth (1992)). 

We assume that U^o is the current time of the particle, 
U,i = ti,o + At^^o is the time after a step is taken, and At^^o 
is the current time step. 

1 . We predict the position and velocity of the particle using 
already known data: 

Xi,pred = x^,o + Ati,oVi,o + -At^o^^,o + ^At^a^,o 



(12) 
(13) 



This operation has a time complexity of 0(1). 
2. We evaluate the acceleration a^ and its derivative (jerk) 
a^ of the particle. To achieve this, we have first to pre- 
dict the positions and velocities of all the other particles 
to time tf. The acceleration and jerk are evaluated on 
either the GRAPE-6 or the CPU using 



ij 



(14) 



(15) 



where 

Yij = Xj^pred ^i^^ved't (1^) 
^ij Vj,pred V^,pred- (1'^) 

If softening e is included in the calculations (for col- 
lisionless simulations), then the denominators of Eqs. 
(14) and (15) are replaced with (rf^- + e^)^^/^^ and 
{rfj + e^)^^^'^\ respectively. 

We also, evaluate numerically higher order derivatives 
of the acceleration using 



-6(a^,o - a^^i) - At^(4a^^o + 2a^^i) 
Atf 

12(£ii,o - ai,i) + 6Atj(ai,o + a^.i) 
At? 



(18) 
(19) 
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Fig. 5. Representation of Hermite 4* order with block time steps. In the figure, a set of 7 particles distributed in 3 time-blocks 
is shown. See text for description. 



3. 



The evaluation of a^^i and a^^i has time complexity 
0{N)^ when it is performed on the CPU and 0(log A/"), 
when it is performed on the GRAPE-6. The other cal- 
culations have complexity 0(1). 

We correct the position and velocity of the particle using 
higher order terms 



At 



At? 



i,pred 



^^^^o + T7da^,o, (20) 



24 



— ^ai,o + ^ai,o. (21) 
6 



120 

At_^ 
~2l 



The complexity of this operation is 0(1). 

Figure 5 shows how particles distributed in 3 time- 
blocks are updated using the H4 scheme. We note that in 
a typical star cluster, the number of time-blocks usually 
used is 10-15. All particles lie initially at the bottom, where 
t = 0. Initially block contains 3 particles, while blocks 1 
and 2 each contain 2 particles. A black dot represents a par- 
ticle after the correction step (a straight blue arrow repre- 
sents the correction step), while a red dot represents a par- 
ticle after a prediction step (a curved red arrow represents 
the prediction step). Force evaluation is necessary for every 



particle, only before a correction step, i.e., force evaluation 
happens for a particle only when it is black. For the force 
evaluation of a particle, a prediction of all other particles to 
the current time of the particle is needed. It is obvious that 
particles in block are considered for predict-evaluate-and- 
correct every dt. Particles of block 1 predict every dt^ but 
subjet to evaluate-and-correct every 2dt. Finally, particles 
of block 2 predict every dt and evaluate-and-correct every 
Adt. The curved brown arrows show how particles are al- 
lowed to move between blocks depending on the time step 
calculated using Eq. (8). A jump to a higher block is al- 
lowed only if it doubles the time step of a particle (particle 
3 jumps from block to block 1). In contrast, a jump to 
any lower block is allowed if necessary (particle 6 jumps 
immediately from block 2 to block 0). 



2.2.3. Binaries and multiples 

During the evolution of the system, close encounters be- 
tween stars may occur. Close binaries also may be formed 
dynamically, after three-body encounters. Those subsys- 
tems require a small time step to be integrated accurately 
in time and this time step may be much shorter than the 
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shorterst time step Atmin used in the N-body code. This 
means that this subsystem has to be specially treated to 
ensure that it evolves in time with the required accuracy 
and speed. The special treatment used for close encoun- 
ters, binaries, and small multiple systems is described in 
the following lines. 

The time-symmetric Hermite 4*^ order integrator 
(Kokubo et al. 1998) is the numerical integrator for bina- 
ries, multiples, and close encounters. This integrator is more 
accurate than the simple H4 integrator with the penalty of 
being a few times slower. Time-symmetry is achieved by ap- 
plying the force evaluation-correction (EC) steps n times, 
after the prediction {P{EC)^) and by appropriate choice of 
the time step. The value n = 3 is used in the code because 
this is sufficient to achieve time-symmetry. Time-symmetry 
is responsible for transforming the H4 scheme into a sym- 
plectic integrator, which is more appropriate for the time 
evolution of small systems such as binaries. Symplectic inte- 
grators (Kinoshita et al. 1990) conserve the characteristics 
of the integrated system (energy, semimajor axis, eccentric- 
ity, period) with high accuracy, paying of course the penalty 
of time. In reality, when using symplectic integrators the en- 
ergy of the integrated system is bounded and changes peri- 
odically with time. In contrast, non-symplectic integrators 
are faster, but the energy of the integrated system increases 
monotonically with time. A brief description of the integra- 
tion method for binaries and multiple sub-systems is given 
at the beginning of Sect. 2.2, between items (i) and (vii). 

The choice of the time step for the {P{EC)^) integra- 
tor is important to preserve time-symmetry. The time step 
at each time t must be such that if the time were to be 
reversed and the particles move backwards, the time step 
choice would remain the same. We followed the logic of the 
algorithm presented in Kaplan et al. (2005) for the choice 
of the time step. According to this algorithm, the time step 
is only allowed to either increase or decrease by a factor of 
2 or remain the same with respect to its last value and an 
increase in the time step is allowed only at even times. In 
this way, the time step remains symmetric and quantized 
as a power of 2. The current time is characterized as either 
even or odd with respect to the last time step taken by the 
system. We apply this algorithm using as a time step crite- 
rion the shortest among the time steps dth^ of the members 
of the BINARY 



the doubling of the previous time step satisfies the criterion 
of being smaller than dt and greater than (it/2, where dt is 
given by Eq. (23). If the criterion is satisfied, then 2Atoid 
is the time step for the next step. If not, then we follow 
the same logic using Atoid as a test time step and follow 
exactly the same steps as if t was odd. 

For this algorithm to work properly, an appropriate 
choice for ryb must be made. This choice needs to be small 
enough for dt not to change by more than a factor of 2 with 
respect to its previous value. In addition, r^b is responsible 
for providing a smooth transition to the time step from the 
N-body code to the binary module of the code, to ensure 
that a particle that has either just joined a binary system 
or just left it does not change its time step by more than a 
factor of 2. This is required because major changes in the 
time step introduce significant errors in the integration. To 
guarantee this behavior, 77b is calculated for every binary 
at the time of its formation. To do this, we determine the 
particle i with the smallest value of |ai|/|ai| and using Eq. 
(22) we determine 77b by requiring that dt = Dtmin/2, which 
is given by Eq. (4) and is the smallest time step of the N- 
body code; it is also the last time step taken by all the 
future members of the binary or the next time step for all 
the particles that leave a binary. In the case of a binary 
with more than 2 members, we use a shared time step for all 
its members. To avoid major changes in the time step crite- 
rion, we use an even smaller value for r^b- The problem here 
is that if a single particle escapes for example from a triple 
system, leaving behind a close binary system, the time step 
of the escaping particle is forced to jump to Dt min even if it 
was much smaller when the particle was still a member of 
the triple. This necessary change in the time step may in- 
troduce errors in the integration. To solve this problem, an 
individual time step must be used for multiple subsystems. 
In this way, all members of the binary would have their 
own time steps and an escaping particle would increase its 
time step by just a factor of 2. 

Finding the binary systems in an N-body system is 
a slow procedure with time complexity 0(A/'^), unless all 
PARTICLE classes are members of a linked list, as in a bi- 
nary search tree, or GRAPE-6 is used for acceleration. Here 
we use the neighbor-module of GRAPE-6 to search more 
rapidly for binaries and close encounters. Each star i with 
mass rrii is associated with a distance 



ai 



(22) 



where 77b is a parameter to control the accuracy, and ai 
and a-i are the acceleration and its derivative for particle i, 
respectively. 

At each time t, we begin by finding whether t is even 
or odd with respect to the last time step, Atoid, taken. If t 
is odd, we maintain the same time step and by using Atoid 
we determine the time step at both the beginning dto and 
the end dti of that time step and find their average 



dt 



dto -\- dti 



(23) 



We constrain the time step used A^oid to be shorter than 
dt and longer than dt/2. If it satisfies this condition, then 
the chosen time step is the one to use in the integration. 
If not, then the time step must be halved. In the case of 
even time, we follow the same procedure, testing whether 



N 



(24) 



where rrig is the mass of the most massive star of the sys- 
tem, N the total number of stars, and Rc\ is given by 
Eq. (5). During the calculation of the force for the star 
i, this distance is sent to the neighbor-module of GRAPE- 
6, which returns the identities of the stars that lie in dis- 
tances r < i^hi from the star i. Those stars are considered 
as "neighbors" of the star i. Then, a search between the 
neighbors is performed to find out whether one or more 
of them is about to have a close encounter or to form a 
close binary system with star i. The criterion for a close 
encounter between the star i and its neighbor j is 



N 

— {mi+mj)Rcu (25) 



(26) 
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where Rij is the separation between stars i and j, m^, and 
rrij their masses, N is the total number of stars, and Vij is 
the relative velocity of the two stars. The choice of the dis- 
tance Rcritij comes from our finding that in an equal-mass 
system, if the distance between two stars is more than Rc\, 
then they can be integrated accurately using the H4 in- 
tegrator with a time step greater than or equal to Dtmin, 
given by Eq. (4), which is the shortest time step allowed 
by the block time-step scheme. If the separation is smaller 
than i^ci, then a smaller time step is needed, so integrating 
the two particles in time with acceptable accuracy, needs 

some special treatment. The factor ^J + ^j) in Eq. 

(25) is used to generalize the criterion for a non-equal-mass 
system. Finally, the distance that is sent to GRAPE-6 
during the calculation of the force for star z, is simply 5 
times the critical distance between this star and the most 
massive star of the system, to ensure that every close en- 
counter between two stars is detected. 

If the above criteria are fulfilled, a binary with the two 
PARTICLE classes as members is created. The two particle 
classes are removed from the N-body system and replaced 
by a new particle that corresponds to a "virtual" star 
represending the center of mass of the removed stars. This 
"virtual" star has mass 



M = m,- 



and is placed at position 



with velocity 



rrii 



rrij 



rrii 



rrin 



(27) 



(28) 



(29) 



and behaves as if it were a normal star in the rest-frame of 
the N-body system. 

A third star k at distance Rk from the center of mass of 
the binary, becomes a member of the binary transforming 
it to a triple system, if one of the criteria below is satisfied: 

1. In all cases, if the distance of the star k from at least 
one of the members of the binary is smaller than 2 / 3i^ci 

2. If the relative motion of the two stars is a simple close 
encounter (e > 1, where e the eccentricity of the sys- 
tem), when 



Rk < \l^{M^mk)Rch 



(30) 



where M is the total mass of the binary and rrik the 
mass of the third star. 
3. If the relative motion of the two stars is a close binary 
(e < 1), then the dimensionless perturbation 7/^, acting 
on the binary from star k 



--Ml)' 

becomes greater than a critical value, i.e., 

7/c > Tcrit- 



(31) 



(32) 



In Eq. (31), M2 = 2m/M, where M is given by Eq. (27), fh 
is the mean mass of the system, i^b is the size of the binary 
defined to be either its semimajor axis for hard binaries 
(i.e., binaries of binding energy greater than the average 
energy of a particle) or the distance R between its members 
in the case of soft binaries (i.e., binaries with binding energy 
lower than the average energy of a particle). 

The same criteria are used when a fourth particle comes 
close enough to a triple system. The critical value of 7 for 
soft binaries is 7crit = 0.125. For hard binaries, we use the 
value 7crit = 0.015625. In an equal-mass system, according 
to Eq. (31), the first value corresponds to a distance Rk = 
2i?b for the third star, while the second to a distance Rk = 
4i?b- This means that in an equal-mass system, a third star 
would become a member of a soft binary when it lies at a 
distance twice the size of the binary. If the binary is hard, 
the critical distance is 4 times the size of the binary. 

Another way of forming multiple sub-systems with num- 
ber of stars n > 4 is when two binaries or multiples ap- 
proach each other. When the distance between their centers 
of mass becomes smaller than three times the sum of their 
sizes, then the two binaries merge forming a larger system 
where all particles interact strongly with each other. The 
size of a binary is defined to be its semimajor axis, for hard 
binaries or the distance between its members, for soft bi- 
naries. 

Termination of a triple, quadruple, or higher order sub- 
system occurs when the following three criteria are met: 

1. When the dimensionless perturbation 7^ acting on the 
particle- member of the sub-system j on the inner binary, 
becomes less than 7crit 



7j < 7crit- 



(33) 



2. The relative velocity of particle j, Vj^ with respect to the 
inner binary is positive, i.e., particle j moves outwards. 

3. The distance of particle j from any of the other par- 
ticles/members of the multiple system, is greater than 
Rc\. 

The so-called "inner binary" is the hardest double sub- 
system inside the multiple (the one with the greatest bind- 
ing energy). 

In some rare cases, two particles-members of a multiple 
system with number of stars n > 4 are removed together 
from the system and form a close binary. 

Termination of a simple binary occurs when the distance 
between its members becomes greater than Rd. Hard bina- 
ries are not terminated, unless another particle becomes a 
member and interacts strongly with their members. The 
outcome of these interactions may be a softer or harder 
binary, a collision, or even a dissolved system. 

Collisions: Another rare case for termination of a binary is 
a collision between its members. A collision between two 
stars occurs when their distance rjn 



'^ij ^ Ri H" Rj •) 



(34) 



where Ri^ Rj the radii of the two stars. 

The stellar radius Ri for a main sequence star is given 

by 



10 



S. Konstantinidis and K. D. Kokkotas: MYRIAD: A new N-body code for simulations of star clusters 




Fig. 6. Perturbers. Stars with 7 > 7pert are considered as perturbers of a binary system. Those are the stars represented as 
empty (red) circles in the above image. The filled (black) circles represent the members of the binary or multiple sub-system. If the 
dimensionless perturbation 7, given by Eq. (31), becomes greater than 7crit, then the star is included into the binary or multiple 
sub-system. More stars could be added forming a multiple system with up to ^ 10 stars. Those systems are usually dispersed 
quickly. On the other hand, a star with 7 becoming smaller than 7pert is no longer a perturber. Those are the stars represented 
with filled (blue) squares. For an equal mass system, we typically set 7pert ^ 10~^ and 7crit = 0.125 if the inner binary is considered 
as "soft", while 7crit = 0.015625 is assumed if it is considered as "hard". 



where is the mass of the star, Rq is the radius, and Mq 
is the mass of the sun. The parameter a depends on the 
mass rrii (Bowers & Deeming 1984) such that 



1 {rrii < O.5M0) 
0.73 {rrii > 0.5Mq) 



(36) 



If the star is a black hole, then Ri corresponds to its 
Schwarchild radius 

Ri = (37) 
where c is the speed of light. 



Perturbers: The internal motion of a binary or multiple sub- 
system is controlled by the internal forces acting between 
its member stars. External perturbations from stars pass- 
ing close to the center of mass of the sub-system are also 
included. A star k is considered as a perturber of a binary 
system if the dimensionless perturbation 7/., acting from 
the star to the binary, is greater than a critical value 



7/c > 7per 



(38) 



A typical value of 10 ~^ is used for 7pert7 but this value is 
usually modified to limit the total number of perturbers 
for a binary system according to the total number of stars. 
During the force calculation, the external perturbations for 
every star-member of a binary or multiple system are calcu- 
lated and added to their internal values. The perturbation 
in the acceleration and first derivative acting on a star 
i, a member of a binary or multiple system, caused by an 
external star j are calculated using 



^i,pert — ^ 



cm, J 7 



(39) 
(40) 



where acmj is the contribution of particle j to the accel- 
eration of the center of mass, acmj is the contribution of 
particle j to the first derivative of the acceleration of the 
center of mass, and a.ij and are the direct contributions 
of particle j to the total acceleration and jerk of particle z, 
respectively. 

Identifying the perturbers of a binary or multiple sub- 
system is done by using the neighbor-module of GRAPE- 
6. During the force calculation for an "virtual" star a 
distance 



Rr. 



1/3 /2mg 



1/3 



d 



(41) 



is sent to the neighbor- module of GRAPE-6. Where rrig is 
the mass of the most massive star in the system, rrii the 
mass of the "virtual" star, i.e., the total mass of the binary 
or multiple subsystem, and d the size of the subsystem i.e. 
the maximum distance between its members. To avoid the 
case of overflow of the neighbor lists of GRAPE-6, each "vir- 
tual" star is sent alone to GRAPE-6, after cleaning the pre- 
vious neighbor lists from GRAPE-6 memory. To avoid over- 
flow and limit the number of perturbers per binary system 
in systems with tens of thousands of stars, Rp. is also mod- 
ified to ensure that ~ 100 neighbors are returned. This is 
achieved by preventing Rp. becoming greater than 2.5 x re, 
where tq is the distance of the sixth nearest neighbor of 
star i, used to find the core radius of the cluster, as we see 
below, and 2.5 x tq is an estimate of the distance of the 
100th nearest neighbor of star i, assuming that the density 
around star i is constant and unchanged from the time of 
the last calculation of tq. This means that the calculation 
of the core radius should be repeated relatively often, to 
estimate the true number density around each star at any 
time. This calculation is usually iterated a few times every 
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crossing time, which is defined to be 

X ^3 1/2 

as described in Appendix A. With the above restrictions, 
the GRAPE-6 memory does not overflow and GRAPE-6 
returns the identities of the "neighbors" of star i, which 
are the perturbers of the corresponding binary system. 
Nevertheless, if an overflow occurs, then the perturbers are 
found after searching among all the stars of the system and 
the result is used to avoid future overflows. 

Figure 6 shows schematically the distribution of per- 
turbers and non-perturbers for a close binary system. 

We note here a couple of definitions used regularly in 
the remainter of the paper. In all simulations we refer to as 
escapers, the star clusters that are assumed to be isolated 
i.e., for which no external galactic field exists. Stars that 
are moving outwards and reach a distance two times the 
initial size of the cluster, are considered as escapers from 
the system. When a star escapes, the energy of the system 
is reinitialized. In addition we refer to as hard binary a 
close binary system whose energy per unit mass h < —1 
or whose semimajor axis a < Rh{mi + mj)/2, where rrii 
and rrij the masses of its members and Rh its half-mass 
radius. In an equal mass system, the limit for a becomes 
a < Rh/N, where N is the number of stars. 



3. Results and performance 

We measure the accuracy and performance of Myriad by 
using it to evolve different systems in time. The most sim- 
ple tests are those where an equal-mass Plummer sphere 
(Plummer 1911) is used as an initial configuration. We 
evolve systems with different numbers of stars and mea- 
sure the accuracy and the speed of Myriad. As a measure 
of the accuracy of the code, we use the relative error in the 
total energy, after evolving the systems for one N-body time 
unit. The speed of the code is measured by the CPU time 
required for the simulations. We also performed simulations 
of systems with equal-mass stars distributed in a Plummer 
sphere and systems with different initial mass functions dis- 
tributed in either a Plummer sphere or according to a King 
density profile (King 1966), up to core collapse. In these 
simulations, we measure the time of core collapse and com- 
pare the results with other codes and with results found in 
the literature. 

In all simulations, the stars are considered as interact- 
ing point particles. The real dimensions of a star are taken 
into account only when it collides with another star. Stellar 
evolution, stellar mass- loss, and binary evolution are also 
not included. If close black- hole binaries are formed, they 
evolve using purely Newtonian dynamics, without any post- 
Newtonian terms in the equations. Finally, the effects of 
the Galactic tidal forces on the star clusters are ignored 
i.e., star clusters are assumed to be isolated. Initial data 
in all simulations were provided by the Star lab software 
environment. 

The host computer used for these tests is an AMD 
Athlon(tm) 64 Processor 3500-h operating at 2.2 GHz with 
2GB of RAM, connected to a 32-chip GRAPE-6 board and 
running Fedore Core 1 GNU/Linux with the 2.4.22 kernel. 



3.1. Accuracy 

There are many different sources of errors in an N-body 
code. Two of them, which can be easily controlled, are the 
choice of integration scheme and its order. Here, as men- 
tioned earlier, the integrator is 4*^ order. If the time step 
were constant and common for all stars, then the errors 
would scale as 0{dt'^). In codes that use block time steps, 
such as Myriad, the parameter that controls the error is the 
accuracy parameter ry is used in the time-step criterion of 
Eq. (8). Another source of errors, which is clear when small 
values of r] are chosen, is the accuracy of the hardware, both 
CPU and GRAPE-6. Myriad uses double precision for all 
calculations on the CPU, while the accuracy of the calcu- 
lation of the accelerations and derivatives is controlled by 
GRAPE-6 (Makino et al. 2003b). 

We performed a series of experiments using Myriad to 
evolve equal-mass Plummer models, with different particle 
numbers N from t = to t = 1 N-body time units. For all 
the integrations, we recorded the cumulative relative energy 
error (AE/E = {Et=i — Et=o)/Eo) and studied its depen- 
dence on the accuracy parameter rj and particle number N. 
For this study, we varied the accuracy parameter rj from 
10""^ to 0.2 and the particle number from 8192(= SK) to 
65536(= 64K), where IK = 1024. The results are shown in 
Figure 7. 

For small values of 77 (77 < 5 x 10~^), the relative energy 
error is almost constant. This is because for those values 
of T] the integration error is small and the total error is 
dominated by the hardware precision. For values of the ac- 
curacy parameter 77 > 5 x 10~^, the error increases with 
r] as expected. When rj > 0.1, the errors are too large, so 
these choices of accuracy parameter are inappropriate for 
simulations. The typical choices for the accuracy parame- 
ter are ( 0.001 < r] < 0.01) and the relative energy error 
between these limits is AE/Eq < 10"^ for TV < 16K and 
AE/Eq < 10~^ for greater values of the particle number 
N. 

The dependence of the error on the particle number N is 
weak and evident only at the lower values of rj. For greater 
values of TV, the relative energy error saturates to slightly 
smaller limits. 

3.2. Performance 

The elapsed time T for a simulation depends on the accu- 
racy parameter r] and particle number N. The elapsed time 
(wall-clock time) as a function of rj and N for the simu- 
lations of Figure 7 is shown in Figure 8. As expected, T 
grows with smaller rj and greater N. We note that smaller 
T] indicates that shorter time steps are used on average. 

To study the speed of Myriad, we performed a set of 
integrations of equal-mass Plummer models with differ- 
ent particle numbers TV from t = to t = 1 N-body 
time units. For all simulations, we recorded the wall-clock 
time as a function of and the results are shown in 
Figure 9, where N is chosen to be between 512 (= 2^) and 
262 144 (= 2^^). The accuracy parameter used for all simu- 
lations was T] = 0.01. In the same figure, the slopes of A/"^, 
N^, and A^log^^g ^ shown for comparison. 

The complete set of calculations inside Myriad has a 
time complexity of 0{N log N). There is only one calcula- 
tion that has a time complexity 0{mN), where 1 < m < N 
is the average number of neighbors per particle returned 
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Fig. 7. Cumulative relative energy error as a function of 77 for a Plummer model with N = 8192 (top left), N = 16384 (top right), 
N = 32768 (bottom left), and = 65536 (bottom right). The integrations ended at t = 1 N-body units. 



by GRAPE-6 or, if GRAPE-6 returns an overflow, calcu- 
lated by the CPU. This is the calculation of the core radius 
of the star cluster. This calculation is repeated every di- 
agnostic time step dt diag- For the experiments of Figure 9, 
this calculation took place only once, before t = 1. As the 
number of stars increases, the calculation of the core ra- 
dius takes a comparable amount of time to the integration 
time. This is why initially the slope in Figure 9 is similar to 
the slope of N\ogN, but as N increases, it tends to more 
closely resemble the slope of N'^. We note that all stars at 
t = begin with the shortest time step Dtmin and then, as 
time passes by, they are slowly distributed in all the avail- 
able time blocks. Because of this, the simulation initially 
operates more slowly. We note that primordial binary sys- 
tems were not included in all runs, while no close binary 
systems were dynamically formed until t = 1 N-body time 
units. At later times, dynamically formed binary systems 
would force the simulation to run slower. 



3.3. Binaries and multiples 

Before testing Myriad in simulations of clusters up to core 
collapse, we investigated the accuracy of the the code in 
handling close binaries and multiples. We recall that the 
algorithm used in the binary module of Myriad is the time- 
symmetric Hermite 4*^ order scheme described in Sect. 
2.2.3 and note that the force calculation is done on the 
CPU. 

Figure 10 shows the cumulative relative error in the 
energy for the first 10 periods of binaries consisting of 
equal mass particles and having diff"erent eccentricities. The 
choice of the accuracy parameter 77b in the time step cal- 
culation and the maximum allowed time step are such that 
for all simulations the initial energy errors are of the same 
order. The number of time steps per period for the bi- 
nary with eccentricity e = 0.19 is 720, for the binary with 
e = 0.51 is 2200, for the binary with e = 0.75 is 2200, and 
for the binary with e = 0.91 is 4300. Figure 11 shows how 
the time step changes according to our time step criterion, 
during the first 2 periods of each of the binaries. As shown 
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Fig. 8. Wall-clock elapsed time as a function of N and r] for the 
simulations of Figure 7. All simulations ended at t = 1 in N-body 
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Fig. 9. Wall-clock elapsed time as a function of the number of 
stars N . All simulations ended at t = 1 in N-body units. The 
accuracy parameter 77 had the same value 0.01 in all simulations. 
The slopes of 0(A^^), 0(A^^), and 0(A^log A^) are also shown for 
comparison. 



in the figure, the time step is quantized as a power of 2. 
Finally, the cumulative energy error for evolving the bina- 
ries for 10^ periods is shown in Figure 12. For a clearer 
visual result, we have printed a fraction of the points in the 
figure. It is obvious that the algorithm used for the dynami- 
cal evolution of close binary systems is time-symmetric and 
preserves the energy since there is no linear growth in the 
energy errors. The energy error oscillates between a lower 
and higher value. The lower error is observed during the 
apastron passage, while the higher value is observed during 
the periastron passage. We would have the same behavior 
for the energy error even if we decrease the number of steps 
per period for the simulations. The error would be greater 
but still have an upper limit and no linear growth with 
time. 

Figure 13 shows the cumulative relative energy error 
in the case of a binary with unequal mass members (mass 
ratio =15) and eccentricity e = 0.89. We followed the sim- 
ulation for 10^ periods. The energy error does not grow 
linearly with time, but instead follows some kind of ran- 
dom walk between some limits. This result shows that the 
binary module of Myriad is capable of evolving binaries 
with high mass ratios and high eccentricities for a long pe- 
riod without significant energy errors. The trajectories of 
the two particles around their common center of mass for 
the full simulation of 10^ periods are shown in Figure 14. 

We finally tested the performance of Myriad in handling 
triple systems by evolving the well known Pythagorean 
three-body system. This system consists of three different 
masses initially located at rest at the corners of a right tri- 
angle whose sides have lengths 3,4, and 5. Each of the three 
masses is located in such a way that its mass is equal to 
the length of the opposite side. Mass m\ = 3 is opposite 
the side of length 3, m2 = 4 is opposite the side of length 



4 and ms = 5 is opposite the side of length 5. The three 
masses interact strongly with each other, and their trajec- 
tories after some time are shown in Figure 16. The final 
result is a hard binary consisting of particles with masses 
7723 = 5 and m2 = 4 and an escaping particle. The evolution 
of the cumulative relative energy error for the simulation 
is shown in Figure 15. We note the two "jumps" in the en- 
ergy error. These "jumps" occur when two of the particles 
come very close to each other and as a result more time 
steps are needed to perform an accurate integration. We 
used a constant accuracy parameter r\\, = 0.0001. The cu- 
mulative relative energy error at t = 100 N-body units is 
2.23776 X 10-^ 



3.3.1. Choice of the parameters 

In Myriad, assuming that a good choice for the accuracy 
parameter r] is 0.01 and adjusting the accuracy parameter 
for binary evolution r^b so that there is a smooth transi- 
tion between the binary module and the Hermite integra- 
tor, there is another free parameter that balances between 
speed and accuracy in collisional simulations up to core 
collapse. This is the number of perturbers for every binary 
or multiple subsystem controlled by the critical dimmen- 
sionless perturbation 7pert- This parameter is discussed in 
Sect. 2.2.3 (Eq. (38)). For the rest of this work, we choose 
7pert = 10"'', which is small enough to ensure that the er- 
rors introduced into the simulations by replacing two or 
more particles that lie close to each other, with their center 
of mass, are not significant. Smaller values of 7pert would 
slightly increase accuracy, but the speed of the code would 
be lowered significantly, since the forces of the perturbers 
on the members of a binary are computed on the CPU and 
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Fig. 10. Cumulative energy error (AE / Eq — (Et — Eq) / Eq) as 
a function of time for the first 10 periods for the simulations 
of equal-mass binary systems with 4 different eccentricities. The 
eccentricities of the binaries are e = 0.19 (red line), e = 0.51 
(green dashed line), e = 0.75 (blue dotted line), and e = 0.91 
(cyan dashed-dotted line). 



Fig. 11. Variation in the time step used for the evolution of the 
binary systems of Figure 10 during 2 periods of each binary. The 
time step parameter ?7b is chosen so that the error in all the bi- 
naries is of the same order. 




Fig. 12. Same as Figure 10, but for the duration of 10^ periods for each binary. 
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Fig. 13. Cumulative energy error (AE / Eq = (^t — Eq)/Eq) as Fig. 14. The trajectory of the two stars in the simulation of 

a function of time for 10^ periods for the simulation of a binary Figure 13. The two stars are orbiting around their common center 

system with eccentricity e = 0.89 and mass ratio mheavy/miight = of mass, which is located at 0(0, 0). 
15. 
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Fig. 15. Cumulative energy error {AE/Eo = {Et - Eo)/Eo) as a Fig- 16- Trajectories of the particles in the simulation of the 
function of time for the Pythagorean three-body system. Pythagorean three-body system. The initial positions of the par- 

ticles are marked with black circles. 
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not on the GRAPE-6. Finally, if the value of 7pert is respon- 
sible for large numbers of perturbers per binary system, this 
value is increased so that no more than 100 perturbers per 
binary are used. 

3.4. Core collapse 

The evolution of a star cluster up to and beyond core col- 
lapse is one of the challenges for any N-body code. During 
core collapse, the density at the center reaches its maximum 
value, while the core radius, defined by the equation 

^^ = V ' ^''^ 

which is described in Appendix A, reaches its minimum 
value. Figures. 23 and 17 clearly show this behavior as a 
star cluster of 1024 equal-mass stars evolves beyond core 
collapse. As the cluster evolves, close binary or multiple sys- 
tems are formed and interact with single stars or other bi- 
nary or multiple systems. Those interactions become more 
frequent and violent as the cluster approaches core collapse, 
because the density of stars around the center of the cluster 
reaches its maximum. All of these encounters have to be re- 
solved with an acceptable accuracy. In addition, some stars 
approach the outer bounds of the system and consequently 
escape the system. 

We tested the behavior of Myriad in evolving star clus- 
ters up to core collapse and compared the results with those 
produced by Starlab or found in the literature. For these 
tests, we used 3 different initial configurations. We initially 
used an equal-mass Plummer model and studied the time 
of core collapse and the evolution of the error in the total 
energy. The same simulation with similar parameters was 
performed using Starlab for comparison. The next test was 
the evolution of a star cluster with an initial mass function 
up to core collapse. The stars were again distributed ini- 
tially according to a Plummer density profile. Finally, we 
performed several simulations of star clusters with an initial 
mass function of a broad range of masses, distributed ini- 
tially according to a King density profile. The result of the 
latter simulations was a measure of the core collapse time 
(tec) as a function of the half-mass relaxation time (trix). 
These results were compared with existing results found in 
the literature. 

3.4.1. Equal-mass Plummer models 

We performed a simulation of an equal-mass Plummer 
model up to core collapse to test the stability, accuracy, 
and speed of Myriad. Figure 17 shows the evolution of the 
Lagrangian radii that corresponds to fractions of the total 
mass of the system and the evolution of the core radius 
of a system, which was initially an equal-mass Plummer 
model with N = 1024 stars. It is obvious that the core ra- 
dius and the smallest Lagrangian radii reach a minimum 
at tec — 340 N-body time units, when core collapse occurs. 
From equation 

which is described in Appendix A, we find that for this 
system the initial half-mass relaxation time is t^i^ = 19.92 



N-body time units, so tec — 17trix, which is close to what is 
expected from Eq. (A. 15) and consistent with the results of 
other codes presented on Table 1 of Freitag & Benz (2001) 
and Anders et al. (2009), where a detailed comparison be- 
tween the two major N-body codes Starlab and NB0DY4 is 
presented. In Figure 18, we compare the results of Myriad 
and Starlab for the evolution of selected Lagrangian radii 
of the same system. The small differences in the Lagrangian 
radii computed by the two codes, may be caused by the 
difference in the escaper-removal criterion and small differ- 
ences in the time-step allocation criteria. Because of these 
differences, the simulation performed with Myriad ended 
with 983 stars remaining in the system, while that com- 
pleted with Starlab contained 973. Figure 19 shows the 
evolution of the instantaneous energy error SE for this sim- 
ulation. In the same figure, the simulation error of Starlab 
is presented for comparison. The error is initially small and 
unaffected by close encounters and binary systems. It grows 
slowly with time, as expected for the H4 integrator. As the 
system approaches core collapse, the interactions at the cen- 
ter become more frequent and violent and the energy error 
is contolled by them. Hard binaries form and their interac- 
tions with single stars or other binaries introduce another 
source of error. After core collapse, the average error be- 
comes smaller, but some peaks, due to strong encounters 
between hard binaries, continue to appear. 

The comparison between Myriad and Starlab is also 
shown in Figure 20, where the respective relative energy er- 
ror ratio Myriad/Starlab is presented. The error of Myriad 
at any time step is smaller by 300 times than the error 
of Starlab, while most of the time the error of Starlab 
is greater. After t ~ 250 and until core collapse, the er- 
ror of Myriad becomes several times smaller than that of 
Starlab. After core collapse, both Starlab and Myriad 
have large instantaneous errors and their ratio scatter from 
10~^ to 300. Figure 21 shows the cumulative relative energy 
error for the simulation as it evolves with time. The same 
error of Starlab again is presented for comparison. It is 
obvious that before core collapse the two codes show iden- 
tical behavior and after core collapse the error of Myriad 
remains smaller than that of Starlab. 

The simulation ended dit t ^ 402 N-body time units. In 
Figure 22, the CPU times of Myriad and Starlab are com- 
pared. For Myriad, the total CPU time and the CPU time 
spent in force calculation and binary evolution are shown. 
Before core collapse, the CPU time of Myriad is constantly 
greater than that of Starlab. This is because of the small 
differences between the two codes in the time step criterion, 
which on average causes Myriad to store more particles in 
small time-step blocks. A search for neighbors also occurs 
every time these particles are updated, which slows down 
the simulation significantly. The CPU time of the Myriad 
simulation until core collapse is about T :^ 1.7 hours, while 
the same time for Starlab is about T 1 hour. After 
core collapse, the CPU time of Starlab becomes longer 
because of the formation of hard binaries. Hard binaries 
lead to an increase in the CPU time of Myriad cause it 
to finish in about T :^ 8 hours. The simulation done with 
Starlab, with similar parameters as those used in Myriad, 
took about 13 hours, spending most of the time in the post- 
collapse part. Finally, in Figure 23 we show the evolution 
of the mass density and the number of stars at the core. As 
expected, the core density increases with time and reaches 
a maximum at the time of core collapse. On the other hand. 
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Fig. 17. Core radius (heavy black line) and Langrangian radii 
containing 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 5%, 
and 3% (blue dashed lines from top to bottom) of the total mass. 
The half mass radius is indicated by a heavier dashed line. Core 
collapse is reached at tec — 340 N-body time units or tec — 17trix 
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Fig. 18. Lagrangian radii containing from top to bottom 90%, 
50%, 10%, and 1% of the total mass. The continuous lines are 
the results produced by Myriad, while the blue dashed lines are 
produced using Star lab. 
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Fig. 19. Error in the energy per quarter of the time unit (6E = 
^(t) — ^(t-0.25)) cLS it evolves in time for the simulation of an 
equal-mass Plummer model of 1024 stars. The simulation ended 
at t 402 N-body time units. The dashed line is the energy error 
of Starlab when it was used for evolving the same system. Before 
core collapse, the energy errors are identical. 



Fig. 20. Respective energy error ratio Myriad/Starlab as a func- 
tion of time. 



18 



S. Konstantinidis and K. D. Kokkotas: MYRIAD: A new N-body code for simulations of star clusters 



time [t^iJ 

10 15 




150 200 



time [t^iJ 

10 




time [N-body] 



time [N-body] 



Fig. 21. Cumulative energy error (AE/Eq = (^t — Eq) / Eq) as 
a function of time for the simulation described in Figure 19. The 
dashed line is the energy error of Starlab, while the heavy black 
line is the energy error of Myriad. 



Fig. 22. CPU time for the simulation as a function of simulation 
time. The continuous line is the total CPU time of Myriad, while 
the heavy dashed line is the total CPU time of Starlab. The 
dashed-dotted line shows the CPU time taken for the force calcu- 
lation and neighbor list creation done on GRAPE-6. The dotted 
line is the CPU time required for the evolution of binary systems. 
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Fig. 23. Core mass-density (in N-body units) (top) and number of stars in the core (bottom). At the time of core 
collapse the core density reaches a maximum value 
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the number of stars at the core drops with time and reaches 
a minimum value at core collapse. 

3.4.2. Plummer models with a Salpeter initial mass function 

A simulation of a Plummer model with N = 16 384 stars 
was performed assuming a Salpeter (Salpeter 1955) ini- 
tial mass function (IMF) and limiting masses of miower = 
O.5M0 and mupper = SM©. In Figure 24 we show the evolu- 
tion of the Lagrangian radii and core radius of the system. 
In Figure 25, we illustrate the evolution of the error in the 
total energy of the system for this simulation. Core col- 
lapse is reached at t 2.2trix, which is much earlier than 
expected for an equal-mass system, given empirically by 
Eq. (A. 15). This is expected because heavier stars tend to 
sink to the core faster and because of that core collapse 
occurs earlier. The time taken for a star of mass m to sink 
to the center from a circular orbit at r Tc is given by 
(Portegies Zwart & McMihan 2003) 



ts = 3.3-trix, (45) 
m 

According to Eq. (45), the most massive stars of the above 
system would sink to the center in a time of about — 
0.73trix- Until t = 2.2trix, most of the massive stars are 
located close to the center and this leads to the observed 
core collapse of the cluster. 

3.4.3. King models with an initial mass function 

We performed simulations of star clusters with initial den- 
sity distribution of a Wq = 6 King model (King 1966) up 
to core collapse. We used a Scalo (Scalo 1986) initial mass 
function (IMF) with lower and upper limits of 0.1 Mq and 
IOOM0, respectively. We varied the number of stars from 
A/' = 6122to N = 12 288 and found the core collapse time 
tec as a fraction of the half mass relaxation time tri^ of 
the systems. The systems chosen are similar to those used 
in Portegies Zwart & McMillan (2003), where the simula- 
tions were performed using the Star lab package including 
stellar mass loss from stellar evolution. The results are pre- 
sented in Table 1. As the core collapse time we consider 
the time of the creation of the first dynamically formed 
hard binary system in agreement with Portegies Zwart & 
McMillan (2003). We performed 18 simulations and our re- 
sult for the core collapse time is tec — (0.17±0.05)trix • The 
result is slightly smaller than Portegies Zwart & McMillan 
(2003), which is tec ^ (0.19 ± 0.08)trix • This difference in 
results may be due to the lack of stellar mass loss and stellar 
evolution in Myriad which, according to Portegies Zwart & 
McMillan (2003), tends to delay core collapse. 

The evolution of Lagrange radii and core radius in one 
of the simulations containing TV = 6122 stars is shown in 
Figure 26. For a clearer visual result, we have applied a cu- 
bic spline interpolation to the core radius. The cumulative 
error in the energy for the same simulation is presented in 
Figure 27. In Figs. 28 and 29, we present the same results 
for one of the simulations containing N = 12 288 stars. The 
cumulative energy error in all simulations remained below 
10-3. 



4. Discussion 

Sources of energy errors There are unavoidable numerical 
errors in all N-body simulations. Here we discuss the source 
of numerical errors in simulations with Myriad and also how 
the total CPU time is divided between the different parts 
of the code. 

The pre-collapse error does not depend on the binary 
evolution since only some short-lived close encounters and 
soft binaries need the special treatment of the binary mod- 
ule of the code. This error depends mainly on the accuracy 
parameter r] as described in Sect. 3.1. The post-collapse er- 
ror is more complex to analyze, because long-living hard 
binaries form and interact with individual particles and 
other binaries. In some cases, multiple subsystems form. 
All of these uncertainties need to be treated very carefully. 
Since the error in the evolution of isolated hard binaries 
and multiples can be very small as described in Sect. 3.3, 
the relatively large errors in the energy after core collapse 
(see Figures 19, 21, 25, 27, and 29) come mainly from the 
interactions of binaries and multiples with the rest of the 
system. The external perturbations in a binary are con- 
trolled by the dimensionless parameter 7pert discussed in 
Sect. 2.2.3. This parameter plays role in the total energy 
error, since if it is large enough, some of the perturbations 
are ignored and systematic errors are introduced into the 
binary evolution and the evolution of stars that lie close 
to the binary. Our choice of 7pert is lO"*", which probably 
needs adjustment to systems containing stars with large 
mass ratios. 

When a binary or multiple sub-system forms accord- 
ing to the rules presented in Sect. 2.2.3, the transforma- 
tion of coordinates from the center of mass of the cluster 
to the local center of mass of the binary or multiple, in- 
troduce another source of energy error. If that binary is 
a soft binary, it is usually required to form close to peri- 
astron passage and is deformed close to apastron passage 
once every period, until it transforms into a hard binary 
or its eccentricity becomes more than 1. This continuous 
transformation of coordinates introduces systematic errors 
into the simulation. Systematic errors are also introduced 
when a third star orbits with high eccentricity around the 
center of mass of a hard binary. According to the rules of 
Myriad, the triple system would have to form and deform 
often depending on the hardness of the hard binary and the 
eccentricity of the third star. In this case, another source 
of error may introduce large errors into the simulation: the 
third star would have to jump from small time steps, when 
it is a member of the triple, to large time steps, when it is 
considered to be an escaper from the triple and evolves as 
a single star in the whole system. 

Finally, it is not clear what the optimal parameters are 
for the construction and destruction of a binary or a multi- 
ple system, because there are many different cases. A mul- 
tiple sub-system for example may have formed around a 
soft binary or hard binary or even a close encounter be- 
tween two stars may occur, and a particle can escape the 
sub- system under different conditions. The parameters used 
in Myriad are those described in Sect. 2.2.3, but they need 
further examination to reduce systematic energy errors. For 
example, it is unclear when an interaction between 2 hard 
binaries can lead to the formation of a multiple sub-system 
with 4 stars and if such a multiple is formed, it is not easy 
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Fig. 24. Core radius (heavy black line) and Lagrangian radii con- 
taining 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 5%, 3%, 
2%, 1%, and 0.1% (blue dashed lines from top to bottom) of the 
total mass. The half-mass radius is indicated by a heavier blue 
dashed line. The initial condition was a Plummer model of 16384 
stars. A Salpeter mass function with mass limits miower = 0.5Mq 
and TTiupper = 5M0 was chosen as the initial mass function. Core 
collapse is reached at tec — 634 N-body time units or tec — 2.2trix- 
The simulation ended at t ^ 660 N-body time units and took ~ 5 
days 
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Fig. 25. Error in the energy per quarter of the time unit (6E = 
^(t) — ^(t-o 25)) as it evolves in time for the simulation of Figure 
24 



Table 1. The core collapse times for different initial models 



Model 


N 


runs 


(trix) [Myr] 


(tec) [Myr] 


tec [trlx] 


6kw6Scalo3 


6122 


11 


5.6 


1.0 


0.18 ±0.05 


12kw6Scalo3 


12288 


7 


13.8 


2.4 


0.17 ±0.06 



Note: The first column shows the name of the model that describes the basic features of the model. The first numbers together 
with k indicate the number of stars in the model. The number following w is the Wo parameter of the initial King density 
distribution. The next set of letters is the name of the initial mass function (IMF) and the last number is the power of 10 of 
the upper-to-lower-mass-ratio of the IMF. For instance 6kw6Scalo3 is the code name of a system of 6A: = 6 x 1024 = 6122 
stars distributed according to a King density profile with parameter Wo = 6 and masses given by a Scalo IMF with 
Supper /'Hiiower = 10^. For all models presented here, the upper mass limit is 100 M© and the lower mass limit is 0.1 Mq. 
The second column shows the number of stars in the simulation. The third indicates the number of runs for the model. 
The fourth column gives the average half mass relaxation time in Myr. The fifth column shows the average core collapse 
time in Myr calculated using Myriad. In the sixth column, the calculated core collapse time in units of trix is presented. 



to determine under the conditions under which it would be 
destroyed forming smaller sub-systems. 

From the above, it is clear that the main source of error 
in the post-collapse part of a simulation, where hard bina- 
ries form in a dense environment, depends on the treatment 
of the interactions between the binaries and their surround- 
ing environment and on the choices we make for the rules 
of creation and dissipation of multiples. The sudden jumps 
in the time step of particles that interact with hard binaries 
are another important source of error. Finally, the dense en- 
vironment itself introduces large errors, because there are 



relatively many stars sharing small time steps and more 
time steps per time unit, introduce large errors in the en- 
ergy. 

CPU time The total CPU time Tcpu taken for a simulation 
with Myriad can be divided into the CPU time needed for 
the different modules of the code 

^CPU — ^BIN ± ^HERM ± ^FORCE ± ^CHECK ± ^DIAG? (^^) 

where Therm is the CPU time for the predictor, correc- 
tor, and finding the time steps of the Hermite 4 scheme. 
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Fig. 26. Core radius (heavy black line) and Lagrangian radii con- 
taining 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 5%, 
and 3% (dashed lines from top to bottom) of the total mass. The 
half mass radius is indicated by a heavier dashed line. The ini- 
tial condition was a King model with Wo = 6 containing 6144 
stars. A Scalo mass function with mass limits miower = O.IM© 
and mupper = lOOM© was chosen as the initial mass function. 
Core collapse is reached at tec — 20.6 N-body time units or 
tec — 0.21trix, which is the time of formation of the first hard 
binary. 



Fig. 27. Cumulative error in the energy {AE/Eq — (Et — Eo)/Eo) 
as it evolves in time for the simulation of Figure 26. 
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Fig. 28. Core radius (heavy black line) and Lagrangian radii con- 
taining 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 5%, 
and 3% (dashed lines from top to bottom) of the total mass. The 
half mass radius is indicated by a heavier dashed line. The ini- 
tial condition was a King model with Wo = 6 containing 12288 
stars. A Scalo mass function with mass limits miower — 0.1 Mq 
and TTiupper = lOOM© was chosen as the initial mass function. 
Core collapse is reached at tec — 37.2 N-body time units or 
tec — 0.21trix, which is the time of formation of the first hard 
binary. 



Fig. 29. Cumulative error in the energy (AE/Eo = (Et — Eo)/Eo) 
as it evolves in time for the simulation of Figure 28. 
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In the simulation of an IK system up to and beyond core 
collapse presented in Sect. 3.4.1, this time is less than 1% 
of the total CPU time. 

The parameter Tforce is the CPU time required by 
GRAPE-6 to calculate the accelerations and identify the 
neighbors of particles. We reduced this CPU time by asking 
GRAPE-6 to search for and return neighbors only for par- 
ticles that have small time steps. We restricted the search 
for neighbors to the particles that belong to the first 3 time 
blocks, without any change in the accuracy of the code. 
For the simulation of Sect. 3.4.1, Tforce — 51%Tcpu. As 
the number of particles in a simulation increases. Therm is 
expected to increase according to the rule N log N. 

The parameter Tbin is the CPU time spent in the evolu- 
tion of binary and multiple systems that form in a cluster. 
It depends on the number of binaries formed, on the time 
step used for their evolution, and on the number of stars 
that are considered as perturbers of a binary system. If pri- 
mordial binaries exist, this CPU time is expected to decel- 
erate the simulation considerably. In the simulation of Sect. 
3.4.1 where no primordial binaries exist, and hard binaries 
formed just before core collapse, Tbin — 42%Tcpu. In the 
same simulation, Tbin before core collapse is insignificant. 

TcHECK is the CPU time needed to check for encounters 
and binaries in the system. We limited this check to parti- 
cles that have small time steps, since the time step given by 
Eq. (8) is an indicator of the dynamics around a particle. 
If at a given time the time step is small, that means that 
there is a probability that the particle is about to undergo 
a close encounter with some other particle, and a search for 
neighbors is required. If the time-step is large, then the par- 
ticle probably lies in a relatively sparce environment and a 
close encounter with other particles is improbable, at least 
at the current time. Tcheck depends on the density of stars 
at the cluster center and for this reason grows for clusters 
that are close to core collapse. In the simulation described 
in Sect. 3.4.1, we found Tcheck — 4%Tcpu. 

Finally, Tdiag is the CPU time required to calculate the 
global parameters of the cluster, such as the core radius, 
remove escaping stars, and write output snapshots to a file. 
The most "expensive" in terms of CPU time is the calcula- 
tion of the core radius of the system, because it evaluates 
the density around each particle. GRAPE-6 is used for this 
calculation, but the process is still relatively slow. For the 
simulation of Sect. 3.4.1, it is Tdiag 

^ 2%TCPU. 

5. Conclusions 

We have developed a new C-h+ N-body code called Myriad. 
The code can simulate the evolution of star clusters with 
excellent accuracy, while its computational speed is satis- 
factory. The accumulated relative error in the total energy 
in the simulation of N = 16 384 stars with a Salpeter IMF 
(^upper = 5M0, miower = O.SMq), cvolvcd up to and be- 
yond core collapse, is smaller than 2.5 x 10~^, while the 
error in the energy per quarter of N-body time unit is 
smaller than 6 x 10~^. The computational time needed for 
this simulation was less than 3 days, which is compara- 
ble to the computational time of most N-body codes, but 
could be improved. The results for the core collapse time of 
star clusters with an IMF or with equal-mass stars are also 
comparable with those found in the literature or produced 
by other codes such as Starlab. In Anders et al. (2009), 
Starlab and NB0DY4 are compared in detail, and Myriad 



successfully reproduces the result for the core collapse time 
of an equal-mass Plummer model consisting of IK parti- 
cles. 

We now summarize the most important innovations in- 
troduced by Myriad: 

1 . The accuracy parameter rj\^ for the evolution of binary 
and multiple sub-systems is adjusted for every one of 
them by the time of its formation, while in other codes it 
is fixed to some value. Variable r^b makes the transition 
of time steps between the H4 scheme and the binary 
module smooth. In the future, 77b may change during 
the binary evolution to retain the accuracy in desirable 
limits. 

2. The algorithms of Kaplan et al. (2005) and Kokubo 
et al. (1998) have been successfully applied to the 
P{ECY scheme for the evolution of binary and mul- 
tiple sub-systems. The error in the energy of these sys- 
tems does not increase linearly with time, but remains 
bounded within some limits in most of the cases. 

3. We have introduced new rules for creating triple systems 
when single stars are approaching binaries. We do not 
use any distance criterion, but the dimensionless pertur- 
bation parameter 7 defined in Eq. (31), which depends 
on the hardness of the binary, the distance of the third 
star and their masses. We use the same parameter to 
remove a star from a multiple sub-system. 

There are 4 important parameters of Myriad that bal- 
ance between speed and accuracy. These are: 

1. The accuracy parameter 77, used for Aarseth's time step 
criterion given in Eq. (8). This is set by the user before 
beginning a simulation with Myriad. A typical value of 
T] is 0.01. 

2. The accuracy parameter 77b for the time step criterion in 
Eq. (22), used for the evolution of binary and multiple 
subsystems. This parameter is adjusted during runtime 
until a smooth transition between the time step of the 
binary and the Hermite integrator is achieved. 

3. The critical value of the perturbation parameter 7crit5 
discussed in Sect. 2.2.3, above which a particle that per- 
turbs a binary system, becomes a member of the same 
binary, forming a triple system. The value of this pa- 
rameter is set in Myriad to be 0.125 for soft binaries 
and 0.015625 for hard binaries. 

4. The critical value of the perturbation parameter 7pert 
above which a particle is considered as a perturber 
of a binary or multiple subsystem (see Sect. 2.2.3: 
Perturbers). A typical value of this parameter that is 
used throughout the present work is 10~^, but it is au- 
tomatically increased if it results in a large number of 
perturbers per binary. The maximum number of per- 
turbers per binary system is set to 100. 

From the results of the present work, we conclude that 
Myriad has the basic structure of a code that could be used 
to model the evolution of star clusters with IMBHs at their 
centers and study the effects of the presence of an IMBH on 
the system and the possibility of its ejection due to encoun- 
ters and collisions with other stellar-mass black holes. To 
achieve this, additional improvements to the code will be 
required so that more physical phenomena are included in 
the simulations. One of the improvements is to include post- 
Newtonian dynamics into the equations for the evolution of 
close binary black holes. Because of the modular structure 
of the code, this improvement could be made by simply 
adding a function under the binary class. This function 
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would calculate the appropriate post-Newtonian terms (up 
to 2.5 PN) for the acceleration and its first derivative, fol- 
lowing the formalisms found in Blanchet (2002), Berentzen 
et al. (2008), and Kupi et al. (2006). Another extension 
that would make the binary class more complete in treat- 
ing binary black holes is the inclusion of the relativistic 
effect of the asymmetric emission of gravitational radia- 
tion and recoil velocity, during the merging of two orbit- 
ing black holes. According to numerical relativity results, 
when two spinning black holes collide, gravitational radia- 
tion could be emitted asymmetrically. This would lead to a 
recoil velocity in the resulting black hole that might be as 
high as 4000A:m/5 (Gonzalez et al. 2007; Campanelli et al. 
2007b,a; Healy et al. 2009; Herrmann et al. 2007), depend- 
ing on the mass ratio of the initial black holes and the 
directions of their spins, but this velocity might be signifi- 
cantly suppressed by the relativistic alignment of the spins 
(Kesden et al. 2010). Tt is not easy, of course to include 
numerical relativity in an N-body code, but semi-analytical 
formulae, coming from fitting between numerical relativity 
results and post-Newtonian theory (Lousto & Zlochower 
2008; Baker et al. 2008; Lousto et al. 2010) to determine 
the direction and magnitude of the recoil velocity. These 
semi-analytical formulas could be easily included in the bi- 
nary module of Myriad, making the code a tool capable of 
reproducing collisions of black holes realistically. Those two 
extensions include effects beyond the Newtonian theory of 
gravitation, and with extensive tests and comparisons to 
theoretical and already published results, will be presented 
in a future paper. After some development. Myriad should 
be able to study the behavior of BHs that come from merg- 
ers in a cluster of stars, the possible changes in the structure 
of the cluster, and the possibility of the escape of IMBHs 
from clusters (Holley-Bockelmann et al. 2007), after colli- 
sions with stellar-mass BHs. 

Additional improvements to the code could be made by 
adding stellar and binary evolution and the infiuence of the 
Galactic tidal field. This would make the simulation more 
realistic, because mass loss from stellar evolution and the 
tidal force of the Galaxy play an important role in the evo- 
lution of the whole cluster. Because of the code structure, 
the addition of stellar evolution may be achieved by simply 
including the appropriate function as part of the particle 
class, while binary evolution could be included as a simple 
function under the BINARY class. 

Finally, improvements could be made to the speed of the 
code, by connecting the particle classes in a binary search 
tree. The search for close encounters and perturbers would 
then become faster and more accurate, since GRAPE-6 
memory overfiow would be avoided. Another change in the 
structure of the code would be to divide the calculations 
between the CPU, GRAPE-6, and the GPU, such that all 
three parts of the hardware cooperate in making the code 
even faster. A version of the code that uses the GPUs in- 
stead of GRAPE-6 to calculate the forces, is under develop- 
ment. Until now, this version has usee the Sapporo library 
(Gaburov et al. 2009) to emulate the GRAPE-6 functions 
on the GPUs and a shared time step for all particles. Results 
of preliminary tests show that, in agreement with the in- 
troductory paper of Sapporo, the speed of the code using 
4 Nvidia Tesla CI 060 GPUs is comparable with the speed 
of the code when GRAPE-6 is used. Parallelization of the 
code, so that it could run on clusters of computers attached 
to accelerating hardware of the GRAPE-6 family or GPUs, 
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is an improvement that needs to be made to be able to 
simulate realistically large star-systems. 

Appendix A: Units and basic equations 

N-body Units: The units adopted in the simulations are the 
usual N-body units^ (Heggie & Mathieu 1986), where G = 
Mt = Ry = 1, and G is the gravitational constant, Mt is 
the total mass of the system, and Ry is its virial radius. In 
these units, the total energy is Et = —1/4. Transformation 
to physical units can be made if the total mass mt and the 
virial radius Vy are known in physical units. If the mass mt 
is given in solar masses [Mq] and in par sees [pc], the 
mass mi of a star i in solar masses 

mi = Mimt [Mq]. (A.l) 

Any distance R in N-body units is transformed to [pc] using 

r = Rry Ipc]. (A.2) 

In addition, there are functions for transforming time and 
velocity from N-body (T, V) to physical units (t, v) 

t = Tr [Myr] (A.3) 

and 

v = VV* [kms-^], (A.4) 

where 

1/2 

V* = 6.557 X 10"^ (^) [kms-^] (A.5) 

and 

/ \ 1/2 

T* = 14.94 (^-^j [Myr]. (A.6) 

The numerical factors 6.557 x 10~^ and 14.94 come from 
expressing 1 x 10-^{GMq/L*) and (i*VG'M0)V2 

in cgs 

units, respectively, where L* is taken to equal 1 pc and G 
and Mq are expressed in cgs units. 

In the following, we follow the definitions of Aarseth 
(2003). 

Ha If- mass radius: The half mass radius rh is defined as the 
radius of the sphere that has its center at the center of mass 
of the star cluster and contains half of the mass of the total 
system. 

Half-mass relaxation time: The equation for the half-mass 
relaxation time is adopted from Aarseth (2003) and Spitzer 
(1987) 

/Arr?\V2 1 

'* = °-''<gs) ma- <^-^) 

where fh is the mean mass and In A is the Coulomb loga- 
rithm (Spitzer 1987) given by 

lnA = ln77V, (A.8) 

where is the number of stars and the factor 7 is usually 
chosen to be 0.4 (Aarseth 2003), but in some cases values 
of around 0.1 may be more appropriate (Giersz & Heggie 
1991). A typical value of In A used in N-body simulations 
is 10. 

^ http : //en . wikipedia . org/wiki/Natural_units#N-body_ 
units 
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Crossing time: The time to cross the cluster is given by 

.3 . 1/2 



tc 



(A.9) 



where m is the raean raass of the systera. In N-body units, 
it is Tcr = 2\/2. Equation (A.9) is derived from 



tcT = 2rv/cr, 



(A.IO) 



where a is the rms velocity dispersion given at equilibrium 

by 



GNm 
2rv 



(A.ll) 



Core radius: The core radius is defined by the equation 

(A.12) 



lEiPl\ri-r,\ 



' C \ v—v O 5 

V Pi 

where pi is the local density around star i given by 



Pi — 4 3 



(A.13) 



and vq is the distance to the sixth nearest neighbor of star 
i, while Td is the position of the density center whose defi- 
nition follows. 



Density center: The density center of the star cluster is de- 
fined to be 

r, = ^^. (A.14) 
Z^i Pi 

Core collapse time: For an equal-mass system, the time for 
core collapse is (Spitzer 1987) 



tec — IStrlx- 



(A.15) 



For a system containing stars with different masses 
distributed according to a King density profile (King 
1966), the core collapse time becomes (Portegies Zwart & 
McMillan 2003) 

tec ^ 0.2trix. (A.16) 

In an N-body simulation, a clearer definition of the ini- 
tial core collapse time is sometimes required since in some 
cases it is unclear when exactly the core shrinks to a min- 
imum size and then starts its expansion. In many publica- 
tions and for more accurate comparison between different 
codes (Anders et al. 2009), the core collapse time is consid- 
ered to be the time of the formation of the first long-living 
hard binary with binding energy greater than lOO/cT, where 
kT is the equivalent to the thermal energy of a gas given 
for a cluster of N stars by 



Ekin = ^NkT 



(A.17) 



and Ekin is the total kinetic energy of the cluster. In other 
publications (Heggie et al. 2006), linear fitting is applied to 
the diagram of Tc versus time and the core collapse time is 
determined by the intersection of two lines, one being the 
fitting line of data before the core reaches its first minimum 



and the other being the fitting line after the core reaches 
its minimum. 

In some other publications, the core collapse time is 
simply found by observing the diagram Tc versus time where 
there is a clear and sharp minimum. 

In this work, the time of core collapse is determined 
both by simply observing the diagram of the evolution of 
core radius with time and by using the formation of the first 
hard binary with binding energy greater than the lOOkT 
criterion. 
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